library(dplyr)
library(ggcorrplot)
library(ggpubr)
library(patchwork)
# library("here")
# library("bookdown")
# library("downloadthis")
library(vegan)
library(plyr)
library(e1071)
library(tidyverse)
library(viridis)
library(GGally)
library(ggrepel)
library(readr)
library(RColorBrewer)
library(oce)
library(plotly)
library(purrr)
library(furrr)
#install_github("jfq3/ggordiplots")
library(ggordiplots)
source("gg_ordisurf_viridis.R")
LOAD OBJECTS IF HAVE RUN ALREADY
geodist_pa<-readRDS(file.path(dataPath,"inputs/Barents_geodist_pa.rds"))
geodist_r6<-readRDS(file.path(dataPath,"inputs/Barents_geodist_r6.rds"))
mds_pa<-readRDS(file.path(dataPath,"inputs/Barents_mds_pa.rds"))
mds_r6<-readRDS(file.path(dataPath,"inputs/Barents_mds_r6.rds"))
#mds_r6_1000<-readRDS(file.path(dataPath,"inputs/Barents_mds_r6_1000rep.rds"))
ep<-0.8 #epsilon (to avoid having to find and run just this line from code blocks where the mds objects were made)
env_sort <- read.csv(file.path(dataPath,"inputs/SPLITS/BARENTS/LDnoSBhi02_env_sort_2022-09-29.csv")) %>% as.data.frame
otu_sort <- read.csv(file.path(dataPath,"inputs/SPLITS/BARENTS/LDnoSBhi02_otu_sort_2022-09-29.csv")) %>% as.data.frame
#must be same length and equal to unique sample number length - the ordered lists are assumed to be directly relatable on a row by row basis
joinedDat<- left_join(env_sort,otu_sort, by=c("SampID"="SampID"))
joinedDat1<-subset(joinedDat, bathy <= -99)
summary(joinedDat1$bathy)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-498.5 -289.4 -251.7 -251.4 -204.5 -100.4
env_sort<-joinedDat1 %>% select(c(2:351))
otu_sort<-joinedDat1 %>% select(c(2,352:491))
env_sort <- env_sort %>% select(-c("BO22_lightbotmean_bdmean",
"BO22_lightbotltmax_bdmean",
"BO22_lightbotltmin_bdmean",
"BO22_lightbotrange_bdmean"))
env_sort_locations<- ggplot(data = env_sort,
aes(x = X.y,
y = Y)) +
theme_classic() +
geom_point(aes(colour = bathy),
size = 2) +
scale_colour_gradient2(low = "red",
mid = "yellow",
high = "green") +
ggtitle("Location of samples in sorted env file")
env_sort_locations
## Removing NAs
otuCompl <- otu_sort[complete.cases(env_sort[, -c(1:which(colnames(env_sort)=="bathy")-1)]), ] # make sure we have the columms
envCompl <- env_sort[complete.cases(env_sort[, -c(1:which(colnames(env_sort)=="bathy")-1)]), ]
## Removing observations with less than 4 OTUs
sel <- rowSums(otuCompl[, -c(1:2)]) >= 4
otuSel <- otuCompl[sel, ]
envSel <- envCompl[sel, ]
## Removing phosphate
#envSel <- envSel %>% select(-phosphate_mean.tif)
dim(otuSel); dim(envSel)
[1] 977 141
[1] 977 346
env_sel_locations<- ggplot(data = envSel,
aes(x = X.y,
y = Y)) +
theme_classic() +
geom_point(aes(colour = bathy),
size = 2) +
scale_colour_gradient2(low = "red",
mid = "yellow",
high = "green") +
ggtitle("Location of samples in sorted env file")
env_sel_locations
envCompl_ccrem<-env_sort%>%filter(!SampID%in%envCompl$SampID)
envCompl_ccrem_locations<- ggplot(data = envCompl_ccrem,
aes(x = X.y,
y = Y)) +
theme_classic() +
geom_point(aes(colour = bathy),
size = 2) +
scale_colour_gradient2(low = "red",
mid = "yellow",
high = "green") +
ggtitle("Location of removed complete case samples resulting in envSel file")
envCompl_ccrem_locations
inv.sel <- rowSums(otuCompl[, -c(1:2)]) < 4
env.invSel <- envCompl[inv.sel, ]
env.invSel_locations<- ggplot(data = env.invSel,
aes(x = X.y,
y = Y)) +
theme_classic() +
geom_point(aes(colour = bathy),
size = 2) +
scale_colour_gradient2(low = "red",
mid = "yellow",
high = "green") +
ggtitle("Location of removed samples with <4 OTUs resulting in envSel file")
env.invSel_locations
# otu_red <- otu1[1:500, ]
# otu <- otu_red
#
# env_red <- env[1:500, ]
# env <- env_red
otu<-otuSel[,-c(1: which(colnames(otuSel)=="Lithodidae")-1)]
env<-envSel[,-c(1: which(colnames(env_sort)=="bathy")-1)]
table(is.na(otu))
FALSE
136780
table(is.na(env))
FALSE
321433
#str(otuSel)
#str(envSel)
# ## Class 1
# otu1 <- subset(otuSel, envSel$SplitRev == 1)
# env1 <- envSel %>% filter(SplitRev == 1)
#
# ## Class 2
# otu2 <- subset(otuSel, envSel$SplitRev == 2)
# env2 <- envSel %>% filter(SplitRev == 2)
#
# ## Class 4
# otu3 <- subset(otuSel, envSel$SplitRev == 3)
# env3 <- envSel %>% filter(SplitRev == 3)
#
# ## Class 4
# otu4 <- subset(otuSel, envSel$SplitRev == 4)
# env4 <- envSel %>% filter(SplitRev == 4)
#
# ## Class 6
# otu6 <- subset(otuSel, envSel$SplitRev == 6)
# env6 <- envSel %>% filter(SplitRev == 6)
#
# ## Class 7
# otu7 <- subset(otuSel, envSel$SplitRev == 7)
# env7 <- envSel %>% filter(SplitRev == 7)
#
# ## Class 8
# otu8 <- subset(otuSel, envSel$SplitRev == 8)
# env8 <- envSel %>% filter(SplitRev == 8)
# ## Selecting
# otu <- otu2
# env <- env2
#env$coords.x1<-as.numeric(env$coords.x1)
#env$coords.x2<-as.numeric(env$coords.x2)
Make function You might have to drop variables that have been imported as character
otu_pa <- decostand(x = otu[, -c(1)],
method = "pa")
# y = ax^w # power transformation formula
dt <- otu[, -c(1:2)] # species data to transform
x_mn <- min(dt[dt > 0])
x_mx <- max(dt)
rng <- 6 # abundance range
w <- log(rng) / (log(x_mx) - log(x_mn))
a <- x_mn^(-w)
# otu_6 <- a * dt[, -c(1:3)]^w
otu_6 <- a * dt^w
range(otu_6)
[1] 0 6
Sys.time()
[1] "2022-09-29 18:25:00 CEST"
dca_pa <- decorana(veg = otu_pa)
Warning: some species were removed because they were missing in the data
print(dca_pa, head=T)
Call:
decorana(veg = otu_pa)
Detrended correspondence analysis with 26 segments.
Rescaling of axes with 4 iterations.
DCA1 DCA2 DCA3 DCA4
Eigenvalues 0.4930 0.3926 0.2994 0.2039
Decorana values 0.8128 0.3817 0.2571 0.2339
Axis lengths 3.6589 4.0737 3.5442 2.8505
## Bray-Curtis
dist_pa <- vegdist(x = otu_pa, method = "bray")
## Geodist
ep <- 0.8 # epsilon
geodist_pa <- isomapdist(dist = dist_pa, epsilon = ep)
saveRDS(geodist_pa,
file = (file.path(dataPath,"inputs/Barents_geodist_pa.rds")))
took 6.4mins with jonatan’s paralell solution (usually ~10mins)
#
# Sys.time()
# d <- 2
# mds_pa <- list()
#
# for (i in 1:100) {
# mds_pa[[i]]<-monoMDS(geodist_pa,
# matrix(c(runif(dim(otu_pa)[1]*d)),
# nrow = dim(otu_pa)[1]),
# k = d,
# model = "global",
# maxit = 2000,
# smin = 1e-7,
# sfgrmin = 1e-7)
# }
# Sys.time()
## monoMDS
# d <- 2
# mds_r6 <- list()
#
# Sys.time()
# for (i in 1:200) {
# mds_r6[[i]]<-monoMDS(geodist_r6,
# matrix(c(runif(dim(otu_6)[1]*d)),
# nrow = dim(otu_6)[1]),
# k = d,
# model = "global",
# maxit = 2000,
# smin = 1e-7,
# sfgrmin = 1e-7)
# }
# Sys.time()
#Jonatan's upgrade for speed - parallel processing of mds reptitions
library(furrr) # Package to run non sequential functions in parallel
#library(purrr)
plan(multisession)
d <- 2
i <-200 # number of reps
List_geodist_pa <- lapply(seq_len(i), function(X) geodist_pa) # makes a list with the input data repeated as many times as reps wanted
start_t<-Sys.time()
Xmds_pa<-furrr::future_map(List_geodist_pa,
function(x) monoMDS(x,
matrix(c(runif(dim(otu_6)[1]*d)),
nrow = dim(otu_6)[1]),
k = d,
model = "global",
maxit = 2000,
smin = 1e-7,
sfgrmin = 1e-7),
.progress = TRUE)
Sys.time() - start_t
can take a few mins
saveRDS(Xmds_pa,
file = file.path(dataPath,"inputs/Barents_mds_barents_pa.rds"))
Make function?
# Loading geodist object
# geodist_nfi <- readRDS(file = "../SplitRev2_geodist_pa_full.rds")
# Loading mds results
# mds <- readRDS(file = "../SplitRev2_mds_pa_full.rds")
## Extracting the stress of each nmds iteration
mds_stress_pa<-unlist(lapply(Xmds_pa, function(v){v[[22]]}))
ordered_pa <-order(mds_stress_pa)
## Best, second best, and worst solution
mds_stress_pa[ordered_pa[1]]
mds_stress_pa[ordered_pa[2]]
mds_stress_pa[ordered_pa[10]]
## Scaling of axes to half change units and varimax rotation by postMDS
mds_best_pa<-postMDS(Xmds_pa[[ordered_pa[1]]],
geodist_pa,
pc = TRUE,
halfchange = TRUE,
threshold = ep) # Is this threshold related to the epsilon above?
mds_best_pa
mds_secbest_pa <- postMDS(Xmds_pa[[ordered_pa[2]]],
geodist_pa,
pc = TRUE,
halfchange = TRUE,
threshold = ep)
mds_secbest_pa
## Procrustes comparisons
procr_pa <- procrustes(mds_best_pa,
mds_secbest_pa,
permutations=999)
protest(mds_best_pa,
mds_secbest_pa,
permutations=999)
plot(procr_pa)
png(file=file.path(dataPath,"outputs/Barents_procrustes_pa.png"), width=1000, height=700)
plot(procr_pa)
dev.off()
# Extracting ordination axis
ax <- 2
axis_pa <- cbind(mds_best_pa$points,
scores(dca_pa,
display = "sites",
origin = TRUE)[, 1:ax])
ggcorr(axis_pa,
method=c("everything","kendall"),
label = TRUE,
label_size = 3,
label_color = "black",
nbreaks = 8,
label_round = 3,
low = "red",
mid = "white",
high = "green")
ggsave(filename = file.path(dataPath,"outputs/Barents_correlationPCAvsNMDS_PA.png"),
device = "png",
dpi=300 )
# Switching direction of NMDS1
# mds_best$points[, 1] <- -mds_best$points[, 1]
dca_r6 <- decorana(veg = otu_6)
Warning: some species were removed because they were missing in the data
print(dca_r6, head=T)
Call:
decorana(veg = otu_6)
Detrended correspondence analysis with 26 segments.
Rescaling of axes with 4 iterations.
DCA1 DCA2 DCA3 DCA4
Eigenvalues 0.4478 0.3736 0.2680 0.2816
Decorana values 0.4613 0.3696 0.2924 0.2686
Axis lengths 4.2345 4.9928 4.7145 3.0913
## Bray-Curtis
dist_r6 <- vegdist(x = otu_6, method = "bray")
## Geodist
ep <- 0.80 # epsilon
geodist_r6 <- isomapdist(dist = dist_r6, epsilon = ep)
saveRDS(geodist_r6,
file = (file.path(dataPath,"inputs/Barents_below100_geodist_r6.rds")))
Jonatan’s parallel solution took 5.4mins, before it took 10mins
## monoMDS
# d <- 2
# mds_r6 <- list()
#
# Sys.time()
# for (i in 1:200) {
# mds_r6[[i]]<-monoMDS(geodist_r6,
# matrix(c(runif(dim(otu_6)[1]*d)),
# nrow = dim(otu_6)[1]),
# k = d,
# model = "global",
# maxit = 2000,
# smin = 1e-7,
# sfgrmin = 1e-7)
# }
# Sys.time()
#Jonatan's upgrade for speed - parallel processing of mds reptitions
#library(furrr) # Package to run non sequential functions in parallel
#library(purrr)
plan(multisession)
d <- 2
i <-200 # number of reps
List_geodist_r6 <- lapply(seq_len(i), function(X) geodist_r6) # makes a list with the input data repeated as many times as reps wanted
start_t<-Sys.time()
Xmds_r6<-furrr::future_map(List_geodist_r6,
function(x) monoMDS(x,
matrix(c(runif(dim(otu_6)[1]*d)),
nrow = dim(otu_6)[1]),
k = d,
model = "global",
maxit = 2000,
smin = 1e-7,
sfgrmin = 1e-7),
.progress = TRUE)
Progress: ────── 100%
Progress: ──────── 100%
Progress: ────────── 100%
Progress: ──────────── 100%
Progress: ────────────── 100%
Progress: ──────────────── 100%
Progress: ────────────────── 100%
Progress: ─────────────────── 100%
Progress: ───────────────────── 100%
Progress: ─────────────────────── 100%
Progress: ───────────────────────── 100%
Progress: ──────────────────────────── 100%
Progress: ────────────────────────────── 100%
Progress: ──────────────────────────────── 100%
Progress: ────────────────────────────────── 100%
Progress: ──────────────────────────────────── 100%
Progress: ────────────────────────────────────── 100%
Progress: ─────────────────────────────────────── 100%
Progress: ───────────────────────────────────────── 100%
Progress: ─────────────────────────────────────────── 100%
Progress: ───────────────────────────────────────────── 100%
Progress: ─────────────────────────────────────────────── 100%
Progress: ───────────────────────────────────────────────── 100%
Progress: ────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────── 100%
Progress: ─────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────── 100%
Progress: ───────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────── 100%
Progress: ───────────────────────────────────────────────────────────── 100%
Progress: ────────────────────────────────────────────────────────────── 100%
Progress: ────────────────────────────────────────────────────────────── 100%
Progress: ─────────────────────────────────────────────────────────────── 100%
Progress: ─────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%
Progress: ──────────────────────────────────────────────────────────────── 100%Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".Warning: UNRELIABLE VALUE: Future (‘<none>’) unexpectedly generated random numbers without specifying argument 'seed'. There is a risk that those random numbers are not statistically sound and the overall results might be invalid. To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe random numbers are produced via the L'Ecuyer-CMRG method. To disable this check, use 'seed=NULL', or set option 'future.rng.onMisuse' to "ignore".
Sys.time() - start_t
Time difference of 1.793691 mins
saveRDS(Xmds_r6,
file = file.path(dataPath,"inputs/Barents_below100_mds_r6.rds"))
# Loading geodist object
# geodist_nfi <- readRDS(file = "../SplitRev2_geodist_nfi.rds")
# Loading mds results
# mds <- readRDS(file = "../SplitRev2_mds.rds")
## Extracting the stress of each nmds iteration
mds_stress_r6<-unlist(lapply(Xmds_r6, function(v){v[[22]]}))
ordered_r6 <-order(mds_stress_r6)
## Best, second best, and worst solution
mds_stress_r6[ordered_r6[1]]
[1] 0.2544858
mds_stress_r6[ordered_r6[2]]
[1] 0.2544858
mds_stress_r6[ordered_r6[100]]
[1] 0.2554883
## Scaling of axes to half change units and varimax rotation by postMDS
mds_best_r6<-postMDS(Xmds_r6[[ordered_r6[1]]],
geodist_r6,
pc = TRUE,
halfchange = TRUE,
threshold = ep) # Is this threshold related to the epsilon above?
mds_best_r6
Call:
monoMDS(dist = x, y = matrix(c(runif(dim(otu_6)[1] * d)), nrow = dim(otu_6)[1]), k = d, model = "global", maxit = 2000, smin = 1e-07, sfgrmin = 1e-07)
Non-metric Multidimensional Scaling
977 points, dissimilarity ‘bray shortest isomap’, call ‘isomapdist(dist = dist_r6, epsilon = ep)’
Dimensions: 2
Stress: 0.2544858
Stress type 1, weak ties
Scores scaled to unit root mean square, rotated to principal components
Stopped after 93 iterations: Scale factor of gradient nearly zero (< sfgrmin)
mds_secbest_r6<-postMDS(Xmds_r6[[ordered_r6[2]]],
geodist_r6,
pc = TRUE,
halfchange = TRUE,
threshold = ep)
mds_secbest_r6
Call:
monoMDS(dist = x, y = matrix(c(runif(dim(otu_6)[1] * d)), nrow = dim(otu_6)[1]), k = d, model = "global", maxit = 2000, smin = 1e-07, sfgrmin = 1e-07)
Non-metric Multidimensional Scaling
977 points, dissimilarity ‘bray shortest isomap’, call ‘isomapdist(dist = dist_r6, epsilon = ep)’
Dimensions: 2
Stress: 0.2544858
Stress type 1, weak ties
Scores scaled to unit root mean square, rotated to principal components
Stopped after 195 iterations: Stress nearly unchanged (ratio > sratmax)
## Procrustes comparisons
procr_r6 <- procrustes(mds_best_r6,
mds_secbest_r6,
permutations=999)
protest(mds_best_r6,
mds_secbest_r6,
permutations=999)
Call:
protest(X = mds_best_r6, Y = mds_secbest_r6, permutations = 999)
Procrustes Sum of Squares (m12 squared): 1.588e-05
Correlation in a symmetric Procrustes rotation: 1
Significance: 0.001
Permutation: free
Number of permutations: 999
plot(procr_r6)
png(file.path(dataPath,"outputs/Barents_below100_procrustes_r6.png"), width=1000, height=700,) #added 1000
plot(procr_r6)
dev.off()
png
2
#### 1000 reps - don’t run if loading saved object Currently commented out as it gained nothing but took extra time. Can be removed in due course.
retain the 200 rep version
# Extracting ordination axis
ax <- 2
axis_r6 <- cbind(mds_best_r6$points,
scores(dca_r6,
display = "sites",
origin = TRUE)[, 1:ax])
ggcorr(axis_r6,
method=c("everything","kendall"),
label = TRUE,
label_size = 3,
label_color = "black",
nbreaks = 8,
label_round = 3,
low = "red",
mid = "white",
high = "green")
ggsave(filename = file.path(dataPath,"outputs/Barents_below100_correlationPCAvsNMDS_r6.png"),
device = "png",
dpi=300 )
Saving 7 x 7 in image
# Switching direction of NMDS1
# mds_best$points[, 1] <- -mds_best$points[, 1]
## Adding scores to data frame
otu_6$gnmds1 <- mds_best_r6$points[, 1]
otu_6$gnmds2 <- mds_best_r6$points[, 2]
otu_6$dca1 <- scores(dca_r6, display = "sites", origin = TRUE)[, 1]
otu_6$dca2 <- scores(dca_r6, display = "sites", origin = TRUE)[, 2]
p_gnmds_r6 <- ggplot(data = otu_6,
aes(x = gnmds1,
y = gnmds2)) +
theme_classic() +
coord_fixed() +
ggtitle("GNMDS",
subtitle = "First run") +
geom_point(colour = "red") +
geom_vline(xintercept = 0,
linetype = 2,
colour = "lightgray") +
geom_hline(yintercept = 0,
linetype = 2,
colour = "lightgray")
p_dca_r6 <- ggplot(data = otu_6,
aes(x = dca1,
y = dca2)) +
theme_classic() +
coord_fixed() +
ggtitle("DCA",
subtitle = "First run") +
geom_point(colour = "red") +
geom_vline(xintercept = 0,
linetype = 2,
colour = "lightgray") +
geom_hline(yintercept = 0,
linetype = 2,
colour = "lightgray")
p_gnmds_r6 + p_dca_r6
NA
NA
ggsave(file.path(dataPath,"outputs/Barents_below100_gnmds_dca_r6.png"),
device = "png",
dpi=300)
Saving 7 x 7 in image
ord <- mds_best_r6
## Axis scores if selected ord is GNMDS
axis <- ord$points %>% as.data.frame
## Axis scores if selected ord is DCA
# axis <- scores(ord,
# display = "sites",
# origin = TRUE)[, 1:ax])
decided to make MLD-bathy vars
env<-env %>%
mutate ("MLDmean_bathy"=MLDmean_Robinson-(bathy*-1),
"MLDmin_bathy"=MLDmin_Robinson-(bathy*-1),
"MLDmax_bathy"=MLDmax_Robinson-(bathy*-1))
env$MLDmean_bathy<-cut(env$MLDmean_bathy,
breaks=c(-2560, -20,20,130),#checked range of values first (min -2554, max 123)
labels=c('belowMLD','onPycno','inMixLayer'))
env$MLDmin_bathy<-cut(env$MLDmin_bathy,
breaks=c(-2560, -20,20,130),#checked range of values first (min -2554, max 123)
labels=c('belowMLD','onPycno','inMixLayer'))
env$MLDmax_bathy<-cut(env$MLDmax_bathy,
breaks=c(-2560, -20,20,130),#checked range of values first (min -2554, max 123)
labels=c('belowMLD','onPycno','inMixLayer'))
env$swDensRob_avs<-swRho(salinity=env$Smean_Robinson,
temperature=env$Tmean_Robinson,
pressure=(env$bathy*-1),
eos="unesco")
env_cont<-env%>% select(-c(landscape,sedclass,gmorph, MLDmean_bathy, MLDmax_bathy, MLDmin_bathy, #categorical
optional, #not a var
MLDmax_Robinson, MLDmean_Robinson, MLDmin_Robinson, #replaced by new vars
MLDsd_Robinson #not meaningful
))
env_cont<-env_cont%>% mutate_if(is.integer,as.numeric)
env_corr <- env_cont # %>% select(-c(SampID))
# env_corr$coords.x1<-as.numeric(env_corr$coords.x1)
# env_corr$coords.x2<-as.numeric(env_corr$coords.x2)
env_corr[(!is.numeric(env_corr)),]
# Vector to hold correlations
cor_ax1 <- NULL
cor_ax2 <- NULL
pv_ax1 <- NULL
pv_ax2 <- NULL
# NMDS1
for( i in seq(length(env_corr))) {
ct.i <- cor.test(axis$MDS1,
env_corr[, i],
method = "kendall")
cor_ax1[i] <- ct.i$estimate
pv_ax1[i] <- ct.i$p.value
}
Warning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zero
# NMDS2
for( i in seq(length(env_corr))) {
ct.i <- cor.test(axis$MDS2,
env_corr[, i],
method = "kendall")
cor_ax2[i] <- ct.i$estimate
pv_ax2[i] <- ct.i$p.value
}
Warning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zeroWarning: the standard deviation is zero
cor_tab <- data.frame(env = names(env_corr),
ord_ax1 = cor_ax1,
pval_ax1 = pv_ax1,
ord_ax2 = cor_ax2,
pval_ax2 = pv_ax2)
cor_tab
write.csv(x = cor_tab,
file = file.path(dataPath,"inputs/Barents_below100_cor-table_r6_200rep_MLD-bathy.csv"),
row.names = FALSE)
cor_a1_sort<-cor_tab%>%
mutate(abs_ord_ax1=abs(ord_ax1),
abs_ord_ax2=abs(ord_ax2)) %>%
arrange(desc(abs_ord_ax1))
cor_a2_sort<-cor_tab%>%
mutate(abs_ord_ax1=abs(ord_ax1),
abs_ord_ax2=abs(ord_ax2)) %>%
arrange(desc(abs_ord_ax2))
dotchart(cor_a1_sort$abs_ord_ax1, main="Absolute (+/-) correlations between envVars and gnmds axis 1")
cor_cut<-0.38 #decide
cor_sel<-subset(cor_a1_sort,abs_ord_ax1>cor_cut)
cor_sel
as.data.frame(cor_sel$env)
env_os <- env[, cor_sel$env]
env_os
str(env_os)
'data.frame': 977 obs. of 20 variables:
$ BO22_dissoxltmax_bdmean: num 329 329 330 328 328 ...
$ BO22_dissoxltmin_bdmean: num 291 291 292 290 290 ...
$ BO22_dissoxmean_bdmean : num 310 310 311 309 309 ...
$ BO22_dissoxrange_bdmean: num 310 310 311 309 309 ...
$ temp_min : num 0.293 0.289 0.188 0.356 0.339 ...
$ BO22_icethickmean_ss : num 1.72e-06 1.77e-06 3.98e-06 5.60e-07 6.01e-07 ...
$ BO22_icethickltmax_ss : num 1.68e-05 1.70e-05 2.95e-05 5.41e-06 5.57e-06 ...
$ BO22_icecoverltmax_ss : num 3.89e-06 3.96e-06 7.49e-06 1.56e-06 1.60e-06 ...
$ BO22_icecovermean_ss : num 5.55e-07 5.82e-07 1.01e-06 0.00 0.00 ...
$ Tmin_Robinson : num 3.77 3.76 3.56 4.1 4.09 ...
$ temp_mean : num 1.59 1.58 1.26 1.86 1.85 ...
$ Y : num 8125334 8125534 8137534 8112534 8112934 ...
$ coords.x2 : num 8125285 8125495 8137580 8112629 8112862 ...
$ BO22_icecoverrange_ss : num 3.97e-05 4.02e-05 7.06e-05 1.62e-05 1.65e-05 ...
$ BO22_icethickrange_ss : num 1.53e-04 1.56e-04 2.74e-04 4.96e-05 5.07e-05 ...
$ Tmean_Robinson : num 4.07 4.06 3.9 4.27 4.26 ...
$ temp_max : num 3.82 3.82 3.92 3.78 3.78 ...
$ Leieschara_sp. : num 0 0 0 0 0 0 0 0 0 0 ...
$ X.y.y : int 485 486 488 490 491 492 493 494 495 496 ...
$ Tmax_Robinson : num 4.54 4.54 4.4 4.62 4.62 ...
ordsrfs <- list(length = ncol(env_os))
for (i in seq(ncol(env_os))) {
os.i <- gg_ordisurf(ord = ord,
env.var = env_os[, i],
pt.size = 1,
# binwidth = 0.05,
var.label = names(env_os)[i],
gen.text.size = 10,
title.text.size = 15,
leg.text.size = 10)
ordsrfs[[i]] <- os.i$plot
}
ordsrfs_plt <- ggarrange(plotlist = ordsrfs,
nrow = 4,
ncol = 5)
ordsrfs_plt
ggexport(ordsrfs_plt,
filename = file.path(dataPath,"outputs/Barents_below100_ordisurfs_top_corr.png"),
width = 2000,
height = 2500)
env_os_m <- env[,c("Tmean_Robinson", #top corr
"salt_max", #top corr
"Smax_Robinson", #comparison to top corr
"swDensRob_avs", #top corr
"BO22_icecoverltmax_ss",#top corr ax2
"BO22_icecovermean_ss",#top corr
"BO22_dissoxmean_bdmean",#top corr
#"BO22_carbonphytoltmin_bdmean",#top corr - no clear gradient in ordisurf
"BO22_ppltmin_ss", #top corr
"X.y", #comparison to Y
"Y", #top corr
"spd_std", #top corr ax2 (blended model)
"CSpdsd_Robinson", #comparison to top corr ax2 (blended model)
"mud", #highest sed var ax1 + corr
"gravel",#highest sed var ax1 - corr
"BO22_silicateltmax_bdmean", #just under top corr ax1
"bathy" #intuitive for comparisons
)]
env_os_m
str(env_os_m)
'data.frame': 1141 obs. of 16 variables:
$ Tmean_Robinson : num 4.07 4.07 4.06 3.9 3.9 ...
$ salt_max : num 35 35 35 35 35 ...
$ Smax_Robinson : num 35.1 35.1 35.1 35.1 35.1 ...
$ swDensRob_avs : num 1029 1029 1029 1029 1029 ...
$ BO22_icecoverltmax_ss : num 3.86e-06 3.89e-06 3.96e-06 7.71e-06 7.49e-06 ...
$ BO22_icecovermean_ss : num 5.82e-07 5.55e-07 5.82e-07 1.06e-06 1.01e-06 ...
$ BO22_dissoxmean_bdmean : num 310 310 310 311 311 ...
$ BO22_ppltmin_ss : num 0 0 0 0 0 0 0 0 0 0 ...
$ X.y : num 1075461 1075261 1075261 1073461 1073261 ...
$ Y : num 8125134 8125334 8125534 8137734 8137534 ...
$ spd_std : num 0.054 0.0541 0.0543 0.0516 0.0516 ...
$ CSpdsd_Robinson : num 0.0222 0.0223 0.0224 0.0163 0.0165 ...
$ mud : num 69.5 60 60 69.5 69.5 69.5 60 60 60 69.5 ...
$ gravel : num 0.5 15 15 0.5 0.5 0.5 15 15 15 0.5 ...
$ BO22_silicateltmax_bdmean: num 6.55 6.55 6.55 6.58 6.58 ...
$ bathy : num -295 -292 -292 -270 -270 ...
ordsrfs_m <- list(length = ncol(env_os_m))
for (i in seq(ncol(env_os_m))) {
os.i_m <- gg_ordisurf(ord = ord,
env.var = env_os_m[, i],
pt.size = 1,
# binwidth = 0.05,
var.label = names(env_os_m)[i],
gen.text.size = 10,
title.text.size = 15,
leg.text.size = 10)
ordsrfs_m[[i]] <- os.i_m$plot
}
ordsrfs_plt_m <- ggarrange(plotlist = ordsrfs_m,
nrow = 4,
ncol = 4)
Warning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite numberWarning: Computation failed in `stat_contour()`:
'from' must be a finite number
ordsrfs_plt_m
ggexport(ordsrfs_plt_m,
filename = file.path(dataPath,"outputs/Barents_ordisurfs_man_sel_domean.png"),
width = 2000,
height = 2000)
## Select if any var should be excluded from envfit (makes less busy to read)
env_os_m_envfit<-env_os_m [,c("Tmean_Robinson", #top corr
"salt_max", #top corr
"Smax_Robinson", #comparison to top corr
"swDensRob_avs", #top corr
"BO22_icecoverltmax_ss",#top corr ax2
#"BO22_icecovermean_ss",#top corr
"BO22_dissoxmean_bdmean",#top corr
#"BO22_dissoxltmin_bdmean",#top corr
#"BO22_carbonphytoltmin_bdmean",#top corr - no clear gradient in ordisurf
"BO22_ppltmin_ss", #top corr
"X.y", #comparison to Y
"Y", #top corr
"spd_std", #top corr ax2 (blended model)
# "CSpdsd_Robinson", #comparison to top corr ax2 (blended model)
"mud", #highest sed var ax1 + corr
"gravel",#highest sed var ax1 - corr
"BO22_silicateltmax_bdmean", #just under top corr ax1
"bathy" #intuitive for comparisons
)]
colnames(env_os_m_envfit)<-c("T", #top corr
"sMx", #top corr
"SmaxR", #comparison to top corr
"swDensR", #top corr
"icecovmax",#top corr ax2
#"icecovav",#top corr
"dissoxav",#top corr
#"dissoxmin",#top corr
#"BO22_carbonphytoltmin_bdmean",#top corr - no clear gradient in ordisurf
"ppltmin", #top corr
"X", #comparison to Y
"Y", #top corr
"spd_std", #top corr ax2 (blended model)
# "CSsd", #comparison to top corr ax2 (blended model)
"mud", #highest sed var ax1 + corr
"gravel",#highest sed var ax1 - corr
"SiLtmax", #just under top corr ax1
"bathy" #intuitive for comparisons
)
## Envfot plot
gg_envfit(ord = ord,
env = env_os_m_envfit,
pt.size = 1)
Error in envfit.default(ord, env, choices = choices, perm = perm) :
missing values in data: consider na.rm = TRUE
ggsave(filename = file.path(dataPath,"outputs/Barents_EnvFit_man_sel_cln_domean.png"),
device = "png",
dpi=300 )
Saving 7 x 7 in image
SampID <- env_sort$SampID
env_vis<-env
env_vis$gnmds1 <- otu_6$gnmds1
env_vis$gnmds2 <- otu_6$gnmds2
env_vis$dca1 <- otu_6$dca1
env_vis$dca2 <- otu_6$dca2
dca_int <- ggplot(data = env_vis,
aes(x = gnmds1,
y = gnmds2)) +
theme_classic() +
coord_fixed() +
ggtitle("Interactive DCA sample ID plot",
subtitle = "First run") +
geom_point(aes(colour = factor(SampID))) +
geom_vline(xintercept = 0,
linetype = 2,
colour = "lightgray") +
geom_hline(yintercept = 0,
linetype = 2,
colour = "lightgray")
ggplotly(dca_int)
Error in `check_aesthetics()`:
! Aesthetics must be either length 1 or the same as the data (1137): colour
Backtrace:
1. plotly::ggplotly(dca_int)
2. plotly:::ggplotly.ggplot(dca_int)
3. plotly::gg2list(...)
4. plotly (local) ggplotly_build(p)
5. plotly (local) by_layer(function(l, d) l$compute_aesthetics(d, plot))
6. plotly (local) f(l = layers[[i]], d = data[[i]])
7. l$compute_aesthetics(d, plot)
8. ggplot2 (local) f(..., self = self)
9. ggplot2:::check_aesthetics(evaled, n)
p_mld <- ggplot(data = env_vis,
aes(x = gnmds1,
y = gnmds2)) +
theme_classic() +
coord_fixed() +
ggtitle("GNMDS coloured by proximity to mixed layer depth",
subtitle = "First run") +
geom_point(aes(colour = MLDmean_bathy)) +
geom_vline(xintercept = 0,
linetype = 2,
colour = "lightgray") +
geom_hline(yintercept = 0,
linetype = 2,
colour = "lightgray")
p_mld
env_vis<- env_vis %>%
mutate(
sedclassName = case_when(
sedclass == "1" ~ "SedCoverR",
sedclass == "5" ~ "Rock",
sedclass == "20" ~ "Mud",
sedclass == "21" ~ "MwBlock",
sedclass == "40" ~ "sMud",
sedclass == "80" ~ "mSand",
sedclass == "100" ~ "Sand",
sedclass == "110" ~ "gMud",
sedclass == "115" ~ "gsMud",
sedclass == "120" ~ "gmSand",
sedclass == "130" ~ "gSand",
sedclass == "150" ~ "MSG",
sedclass == "160" ~ "sGravel",
sedclass == "170" ~ "Gravel",
sedclass == "175" ~ "GravBlock",
sedclass == "185" ~ "SGBmix",
sedclass == "205" ~ "S/MwB",
sedclass == "206" ~ "S/MwG/B",
sedclass == "215" ~ "SGBalt",
sedclass == "300" ~ "HardSed",
sedclass == "500" ~ "Biogenic"
)
)
c25 <- c(
"dodgerblue2", "#E31A1C", # red
"green4",
"#6A3D9A", # purple
"#FF7F00", # orange
"black", "gold1",
"skyblue2", "#FB9A99", # lt pink
"palegreen2",
"#CAB2D6", # lt purple
"#FDBF6F", # lt orange
"gray70", "khaki2",
"maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
"darkturquoise", "green1", "yellow4", "yellow3",
"darkorange4", "brown"
)
p_sed <- ggplot(data = env_vis,
aes(x = gnmds1,
y = gnmds2)) +
theme_classic() +
coord_fixed() +
ggtitle("GNMDS coloured by sediment class",
subtitle = "First run") +
geom_point(aes(colour = factor(sedclassName))) +
scale_colour_manual(values=c25)+
geom_vline(xintercept = 0,
linetype = 2,
colour = "lightgray") +
geom_hline(yintercept = 0,
linetype = 2,
colour = "lightgray")+
guides(colour=guide_legend(ncol=2))
p_sed
env_vis<- env_vis %>%
mutate(
landscapeName = case_when(
landscape == "1" ~ "Strandflat",
landscape == "21" ~ "ContSlope",
landscape == "22" ~ "Canyon",
landscape == "31" ~ "Valley",
landscape == "32" ~ "Fjord",
landscape == "41" ~ "DeepPlain",
landscape == "42" ~ "SlopePlain",
landscape == "43" ~ "ShelfPlain",
landscape == "431" ~ "shallowValley"
)
)
p_land <- ggplot(data = env_vis,
aes(x = gnmds1,
y = gnmds2)) +
theme_classic() +
coord_fixed() +
ggtitle("GNMDS coloured by landscape class",
subtitle = "First run") +
geom_point(aes(colour = factor(landscapeName))) +
geom_vline(xintercept = 0,
linetype = 2,
colour = "lightgray") +
geom_hline(yintercept = 0,
linetype = 2,
colour = "lightgray")+
guides(colour=guide_legend(ncol=2))
p_land
p_gmo <- ggplot(data = env_vis,
aes(x = gnmds1,
y = gnmds2)) +
theme_classic() +
coord_fixed() +
ggtitle("GNMDS coloured by landscape class",
subtitle = "First run") +
geom_point(aes(colour = factor(gmorph))) +
geom_vline(xintercept = 0,
linetype = 2,
colour = "lightgray") +
geom_hline(yintercept = 0,
linetype = 2,
colour = "lightgray")+
guides(colour=guide_legend(ncol=2))
p_gmo
cat_var_plots<-p_mld+p_sed+p_land+p_gmo
ggexport(cat_var_plots,
filename = file.path(dataPath,"outputs/Barents_gnmds_catvar.png"),
width = 1000,
height = 800)
p_gmo <- ggplot(data = env_vis,
aes(x = gnmds1,
y = gnmds2)) +
theme_classic() +
coord_fixed() +
ggtitle("GNMDS coloured by sample",
subtitle = "First run") +
geom_point(aes(colour = factor(SampID))) +
geom_vline(xintercept = 0,
linetype = 2,
colour = "lightgray") +
geom_hline(yintercept = 0,
linetype = 2,
colour = "lightgray")+
guides(colour="none", size="none")
ggplotly(p_gmo)
Error in `check_aesthetics()`:
! Aesthetics must be either length 1 or the same as the data (1137): colour
Backtrace:
1. plotly::ggplotly(p_gmo)
2. plotly:::ggplotly.ggplot(p_gmo)
3. plotly::gg2list(...)
4. plotly (local) ggplotly_build(p)
5. plotly (local) by_layer(function(l, d) l$compute_aesthetics(d, plot))
6. plotly (local) f(l = layers[[i]], d = data[[i]])
7. l$compute_aesthetics(d, plot)
8. ggplot2 (local) f(..., self = self)
9. ggplot2:::check_aesthetics(evaled, n)